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Abstract 
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5-pulses and a suitable phase response curve. 
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I. INTRODUCTION 


Many studies of neural networks and, generally, of eoupled oseillators are based on the assump¬ 
tion that the relevant dynamieal properties ean be reprodueed by restrieting the study to dynamioal 
systems eharaeterized by a single variable: the phase. In spite of its simplieity, this setup is indeed 
able to produee a wealth of nontrivial phenomena, ranging from the synehronization transition [jT]- 
[3l, to self-eonsistent partial-synehronization [|4]-0, and ineluding ehimera states EIHl, to name 
just a few. 

The first sueh model was proposed by Winfree in 1967 to eharaeterize biologieal rhythms 

[lOl. In the weak-eoupling limit, it may reduee to the famous Kuramoto model lUl |2l [TTI . 
that is eurrently mueh used to investigate the synehronization properties of various setups. While 
in the Winfree model the eoupling depends on the absolute value of the oseillator phases, in the 
Kuramoto model it depends sinusoidally on phase-differenees. In faet, the Kuramoto model has 
been generalized to the so-ealled Kuramoto-Daido model lfT2ll . where the eoupling is a generie 
funetion of the phase differenee. 

Independently, yet another elass of oseillators is being investigated: the so-ealled pulse-eoupled 
integrate-and-fire oseillators. Here, a single phase-like variable, deseribing the membrane poten¬ 
tial, inereases linearly until it reaehes a threshold, is thereby reset to some speeifie value, and 
simultaneously triggers the emission of a pulse that is responsible for the mutual eoupling. The 
effeet of the pulse onto the reeeiving oseillator is quantified by the phase response eurve. The 
simplest of sueh models was proposed in the eontext of heart aetivity [fT3ll . but is nowadays quite 
popular in eomputational neuroseienee, where it is widely used to elarify the eolleetive dynamies 
of neural eireuits lfT4l . A similar and mueh used model is the leaky integrate-and-fire (LIF) neuron, 
introdueed by L. Lapieque in 1907, even before physiologieal meehanisms of pulse transmission 
were understood m- There, the membrane potential evolves exponentially rather than linearly in 
time. 

Nowadays, whenever oseillatory phenomena have to be investigated, integrate-and-fire and 
Kuramoto-like models are the most used setups, but it is not elear to what extent the resulting 
phenomenology is typieal of the seleeted model. A prominent example to illustrate the laek of 
a general framework is self-eonsistent partial synehronization (SCPS), a regime where identieal 
oseillators are neither loeked, nor eompletely asynehronous. Kuramoto IH found evidenee of 
SCPS in a network of identieal LIF oseillators in the presenee of noise and delayed (5-pulses. 
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Later, van Vreeswijk observed and analysed this regime in an ensemble of LIF oseillators eoupled 
through smooth pulses and in the absenee of external noise [|5l. SCPS may also arise in the simple 
Kuramoto-Sakaguehi model liTTl (sine eoupling with a phase-shift) but only for a partieular value 
of the phase-shift, when it is marginally stable. The onset of a robust SCPS regime is, however, 
possible in a Kuramoto-Sakaguehi-like setup, under the eondition that the phase-shift parameter 
of the sine funetion depends on the order parameter and the eoupling strength [|T6l . This model 
ean be obtained as a phase approximation of nonlinearly eoupled Stuart-Landau oseillators. 

Another example of differenees among the various setups is emergenee of the irregular eol- 
leetive dynamies in an ensemble of heterogeneous LIFs with delayed (5-pulses ffTTll . The setup 
is superfieially analogous to the Kuramoto ensemble, but ehaotie eolleetive oseillations are not 
possible in the latter model IIT^fT^ . 

In this paper we eompare the various model elasses in the minimal setup of identieal globally 
eoupled oseillators. In order to earry on a meaningful quantitative analysis, three models (A, B, 
and C) are seleeted as follows. Model A is the ensemble of LIF neurons extensively studied in 
Ref. Il20l . By then following [l2T1l . model A is mapped, in the weak-eoupling limit, onto a Winfree- 
type ensemble of oseillators, yielding model B. Finally, model C is obtained as an approximate 
reduetion of model B to a Kuramoto-Daido ensemble. 

Our studies reveal that the seenario emerging from the three models is substantially equivalent 
with a eouple of quantitative diserepaneies whieh eoneern the fully synehronous regime: (i) the 
dependenee of the period on the eoupling strength is different in model A already at the leading 
order; (ii) its stability differs in model C. Finally, the equivalenee between models A,B, and C 
implies that a generie LIF model with pulses of finite width ean be mapped onto a model of 
pulse-eoupled oseillators and ^-like pulses whieh ean be more easily simulated with event-driven 
algorithms. To test this eonjeeture a model of the latter type is introdueed (model D). 

More speeifieally, in seetion II, we introduee the various model elasses, diseuss their mutual 
relationships, and briefly reeall the most eommon asymptotie regimes. Seetion III is devoted to 
a quantitative eomparison of the models A, B, C, and D: in praetiee the analytieally estimated 
stability speetra of the splay and synehronous states, as well as the numerieally obtained features 
of the SCPS are mutually eompared. Seetion IV is devoted to a perturbative analysis of SCPS in the 
Kuramoto-Daido setup. The resulting frequeney of SCPS are found to be in exeellent agreement 
with the numerieal findings. The main results and the open problems are summarized in seetion 
V. Finally, the many teehnieal details related to the stability analysis of the different regimes in the 
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various models are eonfined to five appendiees. 


II. DYNAMICAL REGIMES AND MODEL CLASSES 


As it is well-known, globally eoupled ensembles of identieal oseillators ean exhibit two highly 
symmetrie regimes: (i) a fully synehronized state, where all the oseillators are eharaeterized by the 
same phase at any time and (ii) an asynehronous regime, also ealled splay state, where the phases 
are uniformly distributed. The standard way to quantify the degree of synehronization is via the 
so ealled Kuramoto order parameter 


R 


AT-i 


N 

g27ri(/.j 


( 1 ) 


where N is the ensemble size and j = 1,..., A^, is the proper phase resealed within the unit 
interval. The two above mentioned regimes eorrespond to: (i) R = 1 (fully synehronous regime) 
and (ii) R = 0 (asynehronous regime). 

Besides sueh two extrema, partially synehronized states may be eneountered, whose universal¬ 
ity is less elear. Here below we introduee two major elasses: phase models (whieh inelude the 
Winfree model and the Kuramoto-Daido model) and pulse-eoupled integrate-and-fire oseillators. 


A. Phase-oscillator models 

The dynamies of an autonomous limit-eyele oseillator is often deseribed by a single equation 
for the phase variable. Without loss of generality this variable is introdueed so that it evolves 
aeeording to 

0 = i/ = l/r, (2) 

where u (r) is the frequeney (period) of the oseillation. If the given oseillator weakly interaets 
with its environment (weakness here means that the shape of the limit eyele is not substantially 
affeeted by the perturbation), the phase equation modifies to (see [|2l l22l for details and further 
referenees), 

(j) = u + , (3) 

where xjj is the phase of the foreing, Q is a periodie funetion of both arguments, and g quantifies 
the strength of the foreing or eoupling. Without loss of generality, the eonstant eomponent of Q 
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can be incorporated into frequency v which then becomes ^f-dependent. In many cases Q can be 
represented as 

Q(0,V') = r(0)z(v^), (4) 

where r((/)) is the phase response curve (PRC) and Z{'4)) is the. forcing function. In globally 
coupled oscillators, Z{xIj) can be often expressed as the sum of the contributions of the single ele¬ 
ments, in which case, using the standard normalization g —)■ g/N, one obtains the model structure 
proposed long ago by Winfree to describe biological rhythms [|9l[T0ll. 

fi = u +gT{(j)i)^'^S{(pj) . (5) 

j 

In the weak-coupling limit g v, the interaction, rather than being determined by the absolute 
phases, is determined by phase-differences (see, e.g., [f23ll ). With the help of averaging techniques, 
the model @ can be indeed reduced to the so-called Kuramoto-Daido model llT^ 

^ X] - fj) , (6) 

j 

identified by the single coupling function 

G{i)= + (7) 

Jo 

A brief derivation of this known result [ITTI is sketched in appendix The famous Kuramoto- 
Sakaguchi model ifTTll corresponds to G(^) = sin(—^ -|- /3), where = const. The structure 
of the Kuramoto-Daido model can be further simplified: upon choosing a frame rotating with the 
common frequency u one can get rid of the first term in the right hand side. Moreover, by rescaling 
the time variable, one could remove the explicit dependence on the coupling constant. In order to 
facilitate the comparison with the other models we omit such simplifications. 


B. The Abbott-van Vreeswijk model 

The model consists of N pulse-coupled leaky integrate-and-fire (LIF) units, characterized by 
the scalar variables Ui, i = 1,... iV, all restricted to the unit interval. In the context of neural 
networks, uff) is interpreted as the membrane potential; it evolves according to 

uft) = a-Ui +gE{t) , (8) 
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where a — Ui represents the veloeity field that is assumed to be strietly positive (i.e. a > 1), 
while E{t) is the “mean” field arising from the interaetion with the other oseillators and g is the 
eoupling eonstant. The evolution equation is eomplemented by a resetting rule: onee the potential 
Ui reaehes the threshold value Ui = 1, it is reset to Ui = 0, the neuron fires and a spike is emitted, 
whieh eontributes to the generation of the field E. 

In a globally eoupled system, the field E is the linear superposition of the pulses emitted in 
the past by all neurons. The field dynamies ean be deseribed by an additional, linear differential 
equation, whose Green’s funetion eorresponds to the pulse shape [l24l . In the popular model of 
Abbott and van Vreeswijk [|20l , the neuron firing at f = fo produees the so-ealled a-pulse whose 
shape is 

= (9) 


where t > to, and the eorresponding field equation reads 


E{t) + 2aE{t) + a^E{t) = 

n\t„<t 


( 10 ) 


From now on, the model identified by the Eqs. (8|10) will be referred to as model A. 


C. From the Abbott-van Vreeswijk model to phase models 

For a proper eharaeterization of the splay state with the help of the Kuramoto order parameter 
R, see Eq. Q, it is eonvenient to introduee phase 0 G [0,1) as 

(j) =—ulnlL — u/{a + gu)], (11) 

where u is defined by the implieit formula 

h> =—\nr^[l — 1/{a + gv)] (12) 

and (f){u = 0) = 0, (/>(« = 1) = 1. As shown in [[20l . Eq. ([^ is then transformed to 

^i = V + gV{(t)i)e , (13) 

where e = E{t) — v and 

r(0) = —exp[0/z/] (14) 

a + gv 

is the PRC. In this formulation the field in the asynehronous state is E{t) = v ll20ll and this state 
is eharaeterized by i? = 0. Reeall that cf) is taken modulo one, unless stated otherwise. 
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The model strueture is eompleted by the evolution equation for the field e. Equation ( [T^ now 
beeomes 

2 

e + 2ai + a'^e = ^ ^ [5(f —i/] . (15) 

n\tn<t 

Sinee the sum in the r.h.s. ean be separated into eontributions from N neurons, we write e = 
{1/N) Sj, where Sj = N — v) and is the time of the ith spike (eounted 

baekward starting from time t) emitted by the jth neuron. With this representation we reeognize a 
Winfree-type strueture ([^, with a crueial differenee in that Sj eannot be expressed via the loeal in 
time value of phase, but has its own dynamics. 

In the weak coupling limit, however, the phase of each neuron increases approximately linearly 
in time and the spikes are equispaced [|2T]|, so that t — = t — + (i — l)/i/ = (0^■ + i — l)/z/, 

where (j)j is the phase of the jth oscillator at time t. As a consequence, one can turn the explicit 
time dependence of Sj{t) into a phase dependence, as expected for a Winfree model. By using the 
definition of Ea given in Eq. Q and resumming the corresponding series, one obtains 

S{^) = 

u 


^-a/u 


+ 


1 - e-"A (1 - Q-alvy 


— V . 


(16) 


Eqs. ( 5|14|16 ) define model B. 

Next, we introduce model C: it belongs to the Kuramoto-Daido class and is derived via averag¬ 
ing as an approximation of model B. Eor the forcing function S and the PRC given by Eqs. ( |16|14 ), 
Eq. (|7]) yields the coupling function 


+ 1736 ^/" - 174 , 


(17) 


see appendix 1^ for derivation and Eq. (Bl) for the gn coefficients. The function G is plotted 
in Eig. for some parameter values where SCPS emerges and is stable (please notice that all 
simulations below refer to a = 1.3, while the other parameter values may vary). The coupling 
function G does not reveal any special structure: it has one maximum and one minimum within 
the period. It can be checked, that G(0) = G(l); however, G"(0) ^ G"(l)- The implications of 
such properties are extensively discussed in the next section. 


D, Back to pulse coupled oscillators: a computationally efficient model 

As a corollary of the previous analysis, Winfree-type models characterized by different phase- 
response curves T and different forcing function S', but identical convolution products G (see 
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FIG. 1. The coupling function of model C for different values of the coupling strength: g = 0.02 (solid 
line), g = 0.05 (dotted line), 5 = 0.1 (dashed line), and g = 0.2 (bold line). The other parameters are 
a = 1.3 and a = 6. The corresponding frequency values are u = 0.6986, u = 0.7747, v = 0.7722, and 
= 0.8847. 


Eq. Q) are expected to be equivalent. Among them, it is instructive to consider the model with a 
d-like forcing function and r(0) = (7(0), 


= (18) 

j 

where we have subtracted 1 to ensure a zero average of the forcing function like in the original 
setup. As expected for a Winfree-type model, the argument of the d-function here is the phase. It 
can be transformed into a time-dependent function by substituting 0/i/ —)■ f into the argument of 
the d-function 


0. = z/ - gG{(j)i) + ^(7(0i) ^ 6{t - ti), 


(19) 


3 

where ti is the time when any oscillator is reaching the threshold 0 = 1 This is a standard model 
of d-coupled oscillators with a weakly phase-dependent velocity field. In the following we shall 
refer to it as to model D. 

From a computational point of view it is preferable to change variables, introducing 6i, accord¬ 
ing to 


d0» 


i?(0O 


V - gG{(j)i) 

1^0 


( 20 ) 


so that 9i = uq (with a further adjustment of the PRC that has to be divided by R{(f)i)), while the 
interaction terms would still be easy-to-handle 5-spikes. In fact, since the time derivative between 
the spikes is constant, the simulation of this model does not require a differential-equation solver 
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and can be performed very effieiently. The priee to pay is that 6 is no longer appropriate to 
eharaeterize the splay state, as the eorresponding Kuramoto-order parameter would now differ 
from zero. 

Finally, it is neeessary to eomment about a subtle point: sinee the PRC is negative for 0 = 0, 
the effeet of an ineoming spike on the zth neuron whose phase is just above zero may push it 
baekward below zero. If one interprets 0 as a true phase, this would mean that the i\h neuron is set 
below threshold and thus ready to fire again, a phenomenon that does not happen in the original 
formulation of the model. We should in faet interpret 0 as the membrane potential m in Eq. Q and 
avoid the identifieation of 0 < 0 with 1 + 0 . 


III. MODEL COMPARISON 


A. Splay state 


The splay state and, more preeisely, its stability is the first ground where the three models 
ean be eompared. The stability analysis is performed in the thermodynamie limit N ^ oohy 
introdueing the probability distribution t) of the phases and writing the eontinuity equation 

dP 0(0P) dJ 

dt d(j) d(j)' ^ ^ 

where J is the eorresponding eurrent. 

The three models require different approaehes: for instanee in model A it is neeessary to inelude 
the field dynamies into the analysis, while model C does not require any perturbative expansion. 
In all three eases, however, in the small-^f limit the relevant eigenvalues ean be expressed as (see 
appendix]^ for a detailed aeeount of the ealeulations), 

Hn = ^Tiinu + g5n , ( 22 ) 


where 


(5n = 




0 + - 1 


(23) 


a + gv\ (a + 27rim/)^(l + 27rmi/) ’ 
while the eorresponding eigenveetors are Fourier modes of inereasing frequeney. This result re¬ 
veals a perfeet eorrespondenee among the three models in the weak-eoupling limit. 

In partieular, it is interesting to notiee that the splay state beeomes unstable (along the direetion 
identified by the first Fourier mode) if a exeeeds the eritieal value 


a. 


= —1 + a/ 1 + dTT^ 




(24) 
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FIG. 2. Loss of splay state stability in models A (solid curve), B (dashed curve), C (filled circles), and D 
(dotted curve). The curves for models A,B,D are obtained by numerical study of the correspondent model 
for N = 200 oscillators. The curve C corresponds to the perturbative calculations, see Eq. ( [2^ . The vertical 
line identifies the locus of points numerically analyzed in Ref. Q . 


The loss of stability in model A for = 0.3 was diseovered in Ref. O, where it was shown that it 
eorresponds to the onset of SCPS (see below). Our analysis reveals that this critieal phenomenon 
extends down to the weak eoupling limit and is therefore more general than initially believed. 


In Fig. 1^ we report the bifurcation diagram in the plane {g, a), for a = 1.3. The solid curve, 
obtained by simulating model A for large systems, separates the lower region, where the splay 
state is stable from the upper one, where SCPS is observed. The vertical straight line at = 
0.3 corresponds to the interval of a-values investigated by van Vreeswijk. The dashed curve 


corresponds to the perturbative result (24) as well as to the transition line of model C: it provides 
an excellent approximation even for relatively large g values. 


Quite surprisingly, numerical estimates of the transition line for model B do not reveal appre¬ 
ciable deviations from the perturbative prediction, suggesting that higher order terms are almost 
negligible in the Winfree setup (at least up to ~ 1). The same agreement is observed for the 
d-coupled oscillators in model D. 
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B. Synchronous state 


While considering the synchronous regime, it is instructive to monitor not only stability but 
also the period T of the solution as, contrary to the previous case, it is affected by the coupling 
strength. Let r be the period of the uncoupled system. As follows from Eq. ( [T^ for = 0, 
T = — ln(l — 1/a). 

For model A, by making use of some general formulas derived in [l25ll it is found (see appendix 
B that in the weak-coupling limit the period can be written as T = r + 5T, where 

.,2 




where 


H = 


-{e- - 1) 


z/(l —e 


_ -(a-l)T'i 


{a — 1)(1 — (a — 1)^(1 — 

For models B and C it is instead found that (see again appendix [P]) 


6T = 


gra^ 


u 






(25) 


(26) 


(27) 


These expressions indicate that the agreement between the original FIF setup and Winfree and 
Kuramoto-Daido-type models is not perfect: a difference manifests itself already at the first order 
in g, i.e. 

^ - riT. n 

1.2g. 


B,C 


r ar^ 


Although the discrepancy is not small, it is more on a quantitative than on a qualitative level. 

The stability analysis of the synchronous solution for models A and B (for g 1), (again 
performed in appendix]^ yields the Fyapunov exponent 


A = -g— 


ae 


-(e- - 1) 


1/(1-e-("-^)^) 


(28) 


a [(a — 1)(1 — (a — 1)^(1 — e”""^) 

For a > 1 and g > Q the synchronous solution turns out to be unstable, as it can be appreciated in 

Kg-S 

As for the model C, the stability of its synchronous solution is given by A = gG'{Q), where 
G"(0) is the derivative in the origin (see the appendix]^: here it arises an additional difference. The 
point 0 = 0 is to be identified with 0=1, but the derivative of G(0) in the two points is different: 


in practice, this means that the right derivative differs from the left one; Eq. (28) corresponds to 
the right derivative. The difference between the two derivatives is however somehow irrelevant, as 
it does not affect the sign (at least for our selection of the PRC and pulse shape). 
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FIG. 3. The ratio of the Lyapunov exponent and coupling strength, \/g, versus a and a for the synchronous 


solution of models A, B, see Eq. (281. 


Thus, the perturbative analysis shows that in the limit 1, A the Winfree and Kuramoto- 
Daido models are almost but not perfectly equivalent to the LIF model: the leading correction for 
the period of the synchronous regime differ in models B and C. 


C. Partial synchronization 

Self-consistent partial synchronization has been observed only in a few setups, but the stability 
analysis of the splay state discussed above in this section suggests that this phenomenon might be 
more general than so far believed. In fact, here we show that SCPS arises in all A-D models and it 
can be analyzed perturbatively in the weak-coupling limit. 

A way to spot SCPS is via a nonzero value of the Kuramoto order parameter R. In Fig. it can 
be seen that a transition towards such a regime occurs when the inverse pulse-width a is increased. 
The curves obtained for the four models are rather close to each other, confirming an agreement 
that is expected from the perturbative analysis of the splay state. The more sizable deviations 
concern model A, suggesting that the field dynamics is not entirely negligible. Quite remarkably, 
the outcome of model D is also consistent (see the dotted-dashed curve in Fig. Q, confirming that 
the effect of a smooth pulse shape can be harmlessly transferred to the PRC. 

Let us now identify a signature of SCPS: a difference between the average frequency cu of the 
oscillators (the same for all of them) and the frequency of the mean field 

Q = (0) , where 0 = arg 
where (■) means time average. 

The results are plotted in Fig. (for the same parameter values as in Fig. |^. The two frequen¬ 
cies are reported after subtracting the bare frequency u of the splay state to better appreciate the 
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FIG. 4. Average value of the Kuramoto order parameter for g = 0.1, and a = 1.3: the solid, dashed, 
dotted, and dotted-dashed lines refer to models A, B, C, and D, respectively. 



FIG. 5. (Color online) Average frequencies of the oscillators (top panel) and of the mean field vs a; for 


g = 0.1 and a = 1.3. In both cases, the frequency of the splay state is subtracted, see Eq. (30l, and the 
result is scaled with respect to the coupling strength g. Black circles, red squares, blue pluses, and green 
triangles correspond to models A, B, C, and D respectively. The two solid curves are the outcome of the 


perturbative calculations carried out with model C (see Sec. IV i. 


implication of the transition; i.e. we plot the relative frequeneies 


u = u — u , Cl = Q — u 


(30) 


In the upper panel we see that the oseillator frequeney ui vanishes at the eritieal a-value below 
whieh SCPS disappears. All eurves lie below zero: this means that in SCPS the oseillators are 
slower than in the splay state. In the lower panel, one ean see that the mean field frequeney G is 
smaller than that of the oseillators: this is a typieal signature of SCPS: it means that the oseillators 
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FIG. 6. Evolution of the Kuramoto order parameter in model A for gr = 0.1, a = 5 and a = 1.3. 

“move” faster than their distribution. (Cf. with the results for the nonlinear Kuramoto-Sakaguehi- 
like model in Refs. m, where the oseillators can have any frequency relative to the mean field.) 
At the transition, the value of G coincides with the frequency of the Hopf bifurcation. Once again, 
one can notice a similar kind of agreement among the three models. 

Finally, we plot in Fig. [^the time trace of the Kuramoto order parameter R for the model A and 
an a-value above threshold. There, one can see small periodic oscillations, which are still present 
in model B (data not shown), but completely absent in model C. As explained in the next section, 
this behavior is a consequence of the invariance of the evolution equations under a phase shift. 


IV. PARTIAL SYNCHRONIZATON: A PERTURBATIVE APPROACH 


Within the Kuramoto-Daido setup, the forces depend on phase differences. Accordingly, there 
may exist non-uniform phase distributions that move rigidly in time. They can be viewed as fixed 


points of Eq. (21) in a suitably moving frame. The first example of such a regime was perhaps 
discussed in dH, where the author developed an approximate description of the LIF model in the 
presence of delayed pulses. Here below we show that such states, sometimes referred to as rotating 
waves [12^ , are instances of SCPS. The representation of SCPS as a fixed point allows developing 
a perturbative approach and thereby deriving approximate analytical expressions to be compared 
with the numerics. 


Let us start expressing Eq. ( [2T] ) in a frame that rotates with the (yet unkwnown) frequency G, 
by mapping 0 —)■ 0 — Gi, and then set ^ = 0. By assuming that the velocity field is defined as in 
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Eq. (|^ (for N ^ oo), one obtains 

A 

d(j) 


-n + g #G( 0 -^)P(V^) P( 0 ) 


= 0 , 


(31) 


where is an unknown quantity, to be determined self-eonsistently. Upon integrating the above 
equation, 


P{(j)) = g = eonst, 


(32) 


—U 9 J dijjG^cl) — '0)P('0) 

where the probability flux g is also to be determined. Sinee phases are resealed to the unit interval, 
the flux g eorresponds to the differenee between the average frequeney of the oseillators and that 
of the mean field, 

g = Cj — h = uj — VL. (33) 


In general, there maybe two elasses of solutions of Eq. (32), eharaeterized by r/ = 0 and 7 ^ 0, 
respeetively. In the former ease, the expression in square braekets must vanish. By going in 
Eourier spaee, it ean be easily seen that no sueh probability distribution ean satisfy the eondition 
if all Eourier eomponents 7 ^ 0. On the other hand, whenever = 0, P„ is allowed to be 
different from zero. Sueh distributions are just marginally stable and any arbitrarily small amount 
of noise would smooth them out. The only physieally interesting solutions are those of the seeond 
elass. 

Determining P(0) is not an easy task. Eet us start diseussing the parameter region elose to the 
bifureation point, where deviations from a flat distribution are small. It is eonvenient to rewrite 


Eq. (321 in Eourier spaee. 




27rm0 


' GmPmPn 


27ri(mH-n)^ _ 


= V 


(34) 


m^n 


and to deeompose it into equations for the single eomponents, obtaining 

^Pk 9 ^ ^ GYnPraPk—m g^kO ■ 


(35) 


The simulations reported in Eig. [^suggest that higher order harmonies are inereasingly negligible 
upon approaehing the bifureation. Therefore, we restriet the analysis to the modes k = 1 and k = 2 
(notiee that Pq = 1 for normalization reasons, while Go = 0 by definition, sinee the eonstant term 
of the eoupling funetion is absorbed into the frequeney). Erom the equation for the mode k = 0 
we obtain 

Gl\P,\^ + Gl\p2\^ = {n-g)/2g, (36) 
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FIG. 7. First Fourier components of the phase-distribution for = 0.1 in model C versus the distance from 
the critical point. Circles, squares and diamonds refer to m = 1, 2, and 3, respectively. |Pi| corresponds 
to the Kuramoto order parameter, see Eq. Q- The two solid curves (which scale with exponents 1/2, 1, 
respectively) are the outcome of the analytic calculation. The dashed curve is the outcome of a best fit with 
a slope 3/2. 


where the superseript r means that the real part is being eonsidered. For k = 1 and fc = 2 we find, 


Gi + G2P2 + G1P2 — fi /g , 
G,P^ + {G 2 - (l/g)P2 = 0 , 


(37) 

(38) 


where we have assumed (without loss of generality) that Pi is real (the phase of the solution is 
arbitrary and we ean set the origin as we prefer). 

Let us now imagine that upon variation of the eontrol parameter /i, there exists a transition to 


SCPS for /i = /ic- Sinee P 2 = 0 at the transition, from Eq. (37) it follows that, gGi{jjLc) = (2; 
we eall this speeifie value Gq. Therefore, slightly above the threshold, gGi = Gq + gG'iSfx and 
G = Go + where the prime denotes the derivative with respeet to /r, while dG has to be 


determined. A solution of Eq. (37) is, to the leading order, 

dG — gG'iSfi 


P9 = 


gG2 + Gq 


(39) 


so that now Eq. (381 yields Pf. Next, using that Pi (and thus Pf) is real, we obtain (5G from the 
eondition Im(Pf) = 0, see appendixfor details. As a result, we find that (5G ~ Pf ~ 

A physieally meaningful solution Pi ~ y/SJt, exists for > 0, i.e. above the bifureation point. 


and Eq. (38) implies that P 2 grows linearly and is in general eomplex, meaning that it is shifted 
with respeet to the phase of Pi. Einally, negleeting the term proportional to P 2 in Eq. (36), we 
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determine the last unknown, r]. 


(40) 


Q - r] = 2gGl\Pi\‘^ 

Notice that both and the frequency difference (I — t] depend linearly on the control parameter in 
the vicinity of the bifurcation. 

These perturbative results can be compared with the numerical simulations illustrated in the 
previous section: a plays the role of the control parameter /r. By computing Pi and P 2 for a = 1.3 
and = 0.1 (see Appendix]^, one obtains the data reported in Fig. The two frequencies Cj 
and (1 reveal an excellent agreement with the direct simulation of the three models. Moreover, in 
Fig. 1^ one can see that the theoretical results (see the two upper solid lines) reproduce perfectly 
the behavior of the first two Fourier modes of the phase distribution. 

Away from criticality, many Fourier modes come into play and a perturbative scheme is no 
longer effective. The distribution P(0) can be nevertheless accurately determined by interpreting 
Eq. (|3^ as the fixed point of the recursive relation 


Pn+M) 


_ V _ 

g f G{4> - - (l 


(41) 


This equation shows that g can be determined by imposing the normalization of the r.h.s.. Nu¬ 
merical studies have revealed that generically the recursive procedure either converges to the flat 
distribution P((?)) = 1 or develops nonphysical negative values. We have found that upon tuning 
Cl, one can pass from the former to the latter regime, that are separated by a critical (1 value for 
which the recursive procedure converges to a given shape with some shift. Upon changing the ini¬ 
tial distribution, different phase shifts may be found: the correct solution is the one characterized 
by a zero shift (a true fixed point). Luckily, this objective can be reached by controlling a single 
parameter of the initial distribution: we have found that the most effective one, is the width of the 
distribution itself. Altogether, in spite of the fact that the fixed point is a infinite-dimensional func¬ 
tion, its shape can be determined by tuning two parameters only. The outcome of this procedure is 
shown in Fig. [^for a = 4.7. 


V, SUMMARY AND OPEN PROBLEMS 

In this paper we have performed a quantitative comparison of different model-classes of (phase) 
oscillators. A perturbative analysis of integrate-and-fire oscillators and of the corresponding Win- 
free and Kuramoto-Daido models reveals a substantial equivalence. The stability of the splay state 
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FIG. 8. Probability distribution in the partial synchronization regime for a = 4.7 for the model C. The 
solid curve is the outcome of the recursive procedure discussed in the text, while the pluses refer to direct 
numerical simulations with N = 1000. 

is perfectly reproduced: the whole spectrum of eigenvalues coincides for all of the three models 
up to leading order. As for the synchronous solution, the leading correction to its frequency in 
the Winfree and the Kuramoto-Daido models differs from that found in the LIF model. More¬ 
over, the Kuramoto-Daido model fails to reproduce its stability (left stability differs from the right 
stability as a consequence of a nonanaliticity in the coupling function), although the difference is 
quantitative, but not qualitative. 

The comparison has been extended to the SCPS regime which arises from the splay state 
through a Hopf bifurcation. In this case, a mostly numerical analysis reveals again an excel¬ 
lent agreement among the various models. The largest deviations are observed for the LIF model, 
signaling that the field dynamics is not entirely negligible even in the small coupling limit. 

An important consequence of our comparative studies is the overall evidence that SCPS is 
not specific of integrate-and-fire oscillators, but rather universal, instead. In particular, it is not 
necessary to invoke a dependence on the order parameter, as assumed in [|T^ . 

Furthermore, the mapping of the original LIF dynamics onto a Kuramoto-Daido-type model 
has offered the opportunity to develop a perturbative treatment of SCPS. In fact, in such a setup, 
SCPS corresponds to a uniform rotation of the probability density that can be seen as a fixed point 
in a suitably moving frame and thereby analysed with powerful techniques. 

The actual observation of SCPS in a Kuramoto-Daido setup such as model C opens the question 
of identifying the minimal requisites for its observability. If the coupling function is composed of 
only one harmonics (the Kuramoto-Sakaguchi model), it is known that something similar to SCPS 
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can be observed only in the special case of phase shift equal to 7r/2, where it is anyhow marginally 
stable. In a separate publication we will show that it is sufficient to add a second harmonic to 
observe a stable and robust regime of self-consistent partial synchronization. 


Finally, the good correspondence between model D and the other phase models implies that 
restricting the study to (5-coupled integrate-and-fire oscillators is not a true limitation in so far as 
finite pulse widths can be reduced to such a class by suitably adjusting the phase response curve. 
Such an equivalence has practical advantages, as the former class of models is easier to simulate. 

To what extent the correspondence among the models extends to large coupling strengths is 
also not known: this is another point that is worth to investigate in the future. 
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Appendix A: From Winfree to Kuramoto-Daido 

In the weak-coupling limit, the dynamical changes induced by the coupling occur on long time 
scales compared to the period of the intrinsic oscillations and one can thereby invoke averaging 
techniques. With reference to the model Eq. ([^, it is convenient to expand the coupling term into 
Fourier modes. 




(Al) 


By assuming that only 1:1 resonances matter and retaining the secular terms, i.e. those for which 
m = —n, one obtains 




(A2) 


n 


so that Eq. is obtained since S'_„ = S'„. 
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Appendix B: Derivation of the coupling function of model C 


Using Eq. ([T^ together with Eqs. (7|14) the eonvolution integral ean be written as 


G(^) = y / + B)d^ 


z/2 r ^ z/2 

— e - dxp = —h - —h , 
a a a 


wherey4=(l —e ^,B = e “’’(I—e and a = a + gfz/. Taking into aeeount that-0 + .^ 

shall be understood as taken modulo one, we write 


h = [ 

Jo 




{A'l/j + B)dip 


+ e 


gr (a-l)T / -(a-l)V;r 


(A'l/j + i?]d0 , 




h = . 

Jo Ji-^ 


The further integration is straightforward; it yields Eq. 0. where the eoefficients are given by 
the following expressions, 

h>a^(e^ — 1 ) 


9i = 


92 = 


93 = 


(a + gv)(a — l)(e"'^ — 1) ’ 


1 


+ 


z/ 


1 - e-“^ a - 1 ’ 




(a + gu/Za — 1)2 ’ 


gi = 


iy\e^ - 1 ) 

a + gu 


(Bl) 


Appendix C: Linear stability of the splay state 


1. Model A 


The weak-coupling limit of the splay state in this setup has been first studied in 
reeently extended to a broader elass of pulse-eoupled integrate-and-fire systems in 


from Eq. (21) with the flux 


J(0,f) = [u + ge^/’'e{t)] P((f)A) 


and more 
. We start 


(Cl) 
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and with the boundary condition J(0, t) = J(l, t). The evolution equation for the field is 


e{t) + 2ae{t) + a^e{t) = J(l, t) — Eq) . 


(C2) 


The splay state corresponds to Pq = I, ^ = 0, and Jq = Eq = u. 

Upon introducing the perturbation j{(p,t) to the steady flux Jq = u, i.e. writing J{(p,t) = 
V + j{(j), t), the evolution equations ( 21|C1|C2 ) can be linearised, yielding 

^ 


dt a ^ dt d(j) ' 
e{t) + 2ae{t) + a^e{t) = . 


(C3) 

(C4) 


Using the standard Ansatz j{(p,t) = j/(0) exp(/xt) and e{t) = e:/exp(/xt), and imposing the 
boundary condition j/(0) = one obtains the eigenvalue equation 

., 2 ,. ^1 


- 1) (/u + af = . 


(C5) 


We now investigate the weak coupling limit g 1. For g = 0, two eigenvalues are obtained 
by solving (/i + a)^ = 0, i.e. /i = —a is a double degenerate solution. Besides, the spectrum 
consists of an infinite set of purely imaginary eigenvalues, /i = 2Tiinu, n ^ 0, which are most 
important for determination of stability. In the small g limit one can assume = 2'Kinv + g5n- 


Upon replacing in Eq. ( |C5[ ), we obtain 

5n{2,T:inv + aY = 


27rm«V2 f\^^il/u+2nin)^ _ 


a 


>0 


Computing the integral, one obtains the final Eq. (23). 


2. Model B 

Here, we refer to model ([^. In the thermodynamic limit, the sum over all oscillators transforms 
into an integral, and the expression for the probability flux takes the form 

J(0, t) = [u + gnc(>)Sp{t)] P{4>, t) , (C6) 

where 

Sp{t)= [ d'ipS{'il;)P{i;,t) , (C7) 

Jo 

while the boundary condition reads 

[z/ + gr{l)Sp{t)]P{l, t) = [u + gr{0)Sp]P{0, t) . (C8) 


21 







At variance with the previous ease, the stability ean be assessed by just linearizing the above 
equation, without the need of ineluding the field dynamies. The problem ean be formally solved 
for arbitrary eoupling strength 

Starting from Eqs. ( C6|C7 ) with the boundary eondition (C8), we set = 1 + 

where p(0, t) represents a perturbation around the homogeneous solution. The linearized equation 
writes 

f)n fir) 

(C9) 




where the prime denotes derivation with respeet to 0 and Sp is defined analogously io Sp, see 


Eq. (C7); notiee also that S'p = 0 in the splay state. The boundary eondition beeomes 


v[p{l,t) -p{0,t)] = -gSpAT , 


where AT = r(l) — r(0). 

Next, we introduee the usual Ansatz, p(0, t) = p(0)e^^ obtaining 

= -pp - ^r'(0)^p, 


(CIO) 


where Sp is defined analogously to Sp, see Eq. (C7). By assuming that p(0) = po(0) exp(—/r^/ u), 
we find that 


Po{4>) =-^S.iM) + C , 


where 


IM = / . 

Jo 

The integration eonstant ean be determined from the boundary eondition 

g e->/%(l)-Ar 

- 1 ^ ■ 


(Cll) 


(C12) 


u 


As a result. 


p(0) = 


- AT 

Q-ir/u _ 1 


- W) 


Sn 


(C13) 


(C14) 


The eigenvalue equation is finally obtained by multiplying p(0) by 5(0) and integrating over 0 to 


obtain Spi 
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(C15) 


Q-^r/y _ 1 

where (■) 5 denotes the integral over the dummy variable 0 after having been multiplied by 5(0). 


22 











In the weak coupling limit, the second addendum in the l.h.s. of the above equation can be 
neglected, while the first one can be properly handled by assuming = 2Timv + g5n in the 
numerator (and /i„ = 27rinz/ everywhere else). As a result, the eigenvalue equation simplifies to 


= 


f' - Ar 




(C16) 


since Eq. (C12) reduces to the Fourier transform of E', while (e reduces to the conjugate 

of the transform of S. From Eq. ([T4|), it follows that 


/■l y — 1 

r„= / d 0 r( 0 )e 2 ™'^ = 


a l/i/ + 27rm ’ 


and, accordingly. 


so that 


By further noticing that 


f' = 1 ~ ^ 


r; - AF = 


a l/z/ + 27rin ’ 

— 1 27rmz/ 
a l/z/ + 27rm 






{a — 2nini>Y 


(C17) 


(C18) 


(C19) 


(C20) 


we finally obtain Eq. (23). 


3. Model C 


In the thermodynamic limit, Eq. ([^ can be written as 

^ = u + gj d'0G(0 - '0)P(^) , 

or, using the Fourier representation, as 


= u 


gJ^G^Pn 


27rm0 


Accordingly, the continuity equation becomes 


dP 


A 

d(j) 


u 


+ gJ^GnPne-^^'^^A P{cl),t) 

n / 


(C21) 


(C22) 


(C23) 


We now linearize Eq. (C23) around the splay solution Po(0) = 1. by assuming P(0, f) = 1 + 
p(0, t). Since the mode amplitudes of the equilibrium solution for n ^ 0 are all equal to zero. 


dp 2^i^ ^ 

nf^0 


dt 


d(j) 


(C24) 
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At variance with the previous setups, one ean easily solve the eontinuity equation by just going in 
Fourier spaee, as this ehange of variables diagonalizes the evolution equation for any parameter 
value 


dpn 

dt 


2mn 


V + gGn 


Pn ■ 


(C25) 


By reealling that = 27rini/ + g6n. 


6n = 27r'mGn = 27rinr„S'* , 


(C26) 


whieh, in the ease of the LIF model, eoineides with Eq. (|2^. 


Appends D: Linear stability of the synchronous state 


1. Model A 

From Eq. (52) in Il25l . the period T for the an ensemble of EIE oseillators with a-pulses is 
determined by the implieit eondition 


a(l-e '^)+g 


q-T _Q-aT 


where 


Q = 


a — 1 
c? j {a — 1) 


(1/ + Q)-Te-“'Q 


= 1 


^-aT 


V = 


(1 


^-aT\2 


(Dl) 


(D2) 


Eor g = 0, the period is equal to r = — ln(l — 1/a) = 1/z/, of. Eq. (12); let us denote with Qo, Vq 
the oorresponding values of Q and V. In the small g limit, we ean assume T = r + 5T, where 6T 


is small, and expand the first term in Eq. ( |D1[ ) (the seeond term is already of order g), obtaining 

(D3) 


ST=-- 


l_e-(a-l)r 


a 


(K, + go)-re-(“-i)"gc 


By replaeing the expressions for Vq and go, we finally obtain Eqs. ( 25|26 
The stability of the limit oyole is determined by the exponent 


A = - In 


a + gV 


- 1 


T a — 1 + gV 
By now expanding for ^ 1, we obtain, up to the first order, 

gVo 

Ta{a — 1) 


r 


(D4) 


(D5) 


With the help of Eq. (25) and recalling that e’" — 1 = l/(a — 1), one obtains the final expression 
for the Eyapunov exponent that is reported in Eq. 
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2. Model B 


Here, we determine the period and determine the stability of the fully synehronous solution 
of the model Q. In the weak eoupling limit, the period ean be estimated through a perturbative 
ealeulation, by setting cj) = ut + (]{t) in Eq. (j^ and retaining the leading order, 

/3 = gT{i't)S{i't). (D6) 

The period T ean be then obtained by solving the above equation and imposing 

z/T + /?(r) = l, (D7) 


so that 


6T = 


I3{t) ean be determined by integrating Eq. (D 61 that ean be written as, 

■'2 


f3 = 


ga 


+ 




u 


(D8) 


(D9) 


where we have used that d = a in the weak eoupling limit. By replaeing the integral of this 


equation into Eq. (D 81 , one obtains the expression reported in Eq. 
As for the stability, the tangent spaee evolution writes 


= gT'{(j)i){S)5(j)i + gT{(pi)^ ^ . 


df 


(DIO) 


If all the oseillators are synehronized, we ean drop the index dependenee in the phase spaee dy- 
namies. 


^ ^ gr'(4.)Sm4>i + ' <D11) 

j 

The stability ean be assessed by introdueing the variables 9i = 5(f)i — 5(f)i with ^ > 2 1123 (the 
sum of all 5(l)i gives a missing equation whieh is known to yield the zero exponent and we thereby 
avoid eonsidering it). 


e, = gr'{<p)s{<p)ei . 


(D12) 


Sinee r(0) is diseontinuous for 0=1, its derivative has a delta eontribution that has to be properly 
ineluded in the eomputation of the Eloquet exponent. The final result is 


T 0(0) T + (71(1)5(1) ■ 


(D13) 
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In the limit g 1, taking into account that S'(O) = S'(l) the above equation reduces to 

1 


A = - In 


r 0 ( 0 ) 


+ G[r(o)-r(i)]^(o). 


(D14) 


One can then determine 0(r) by integrating Eq. (D12) with the same philosophy as for Eq. (D 6 ). 


As a result the same expression as (28 1 is obtained for A. 


3. Model C 

The determination of the period is pretty straightforward: it can be obtained by setting the 
argument of the interaction function G equal to zero 


— — u + gG{0) , 

so that, for the EIE oscillators, 

5T = -gG{0)T‘^ = g{gig2 + 5'3 - 5'4)r^ 


(D15) 


(D16) 


Upon replacing the expressions for gi, g 2 , g^ and (74 reported in appendix]^ one can verify that 
the above equation coincides with Eq. ( |27] ). 

Next, we linearize the equations of motion, obtaining 

d6(j)i 


dt 


gG\0)S4>i - ^G'(O) ^ . 


(D17) 


The stability can be determined by again introducing the variables 0* = 50* — 50i, which satisfy 
the following equation 

0, = ^G'(O)0, , (D18) 

so that the stability is controlled by the sign of G'( 0 ). 

Appends E: Computation of the first Fourier mode of the prohahility distribution in SCPS state 


Substituting Eq. (38) into Eq. (391, we obtain 


~ 2 _ Uo — gG2 — gG[6ix 
-'1 ~ 


^0 + 0G2 


Q(] 


(El) 


Condition Im(Pf) = 0 yields 

6Q = 


n‘^\n |2 _ 02 

{G[Y + {G[y ^ ' ° 6fi = MSfi 

2gnoGi 


(E2) 
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As a result, from Eq. (El), it follows that Pf is proportional to d/j.. 


2 ^0 ~ 9 G 2 M — gG[ 


Pt = 




(E3) 


^0 + gG2 fio 

To complete the computation we have to find G 2 , G[ at the bifurcation point. With the reference 
to the Abbot - van Vreeswijk model, the coupling function is given by Eq. ( [T7| ) and the role of 
the order parameter is played by the inverse pulse width a; the bifurcation value Oc is given by 
Eq. ( |2^ . Computing the first two Eourier modes of G, we find: 

1^1,2 — iPl,2 


Gi ,2 = C- 


-Di ,2 


(E4) 


where 


and 


G = 




a + gu 


1 ) 


Einally, 


An = al — {27inh'Y{l+2ac) , 

Bn = 2Tmv \a^^ + 2ac —{2TmvY~\ , 
Dn = [ttc + (27rnz/)^]^[l + (27rnz/)^] . 


{G',r 


2G fAi 2acAi 

Di \ etc + ct^ 

4 G 7 rz/(l + etc) 

A ■ 
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